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The spectral expansion of the Green's tensor for a planar multilayered structure allows us to semi analytically 
obtain the angular spectrum representation of the field scattered by an arbitrary dielectric perturbation present 
in the structure. In this paper we present a method to find the expansion coefficients of the scattered field, 
given that the electric field inside the perturbation is available. The method uses a complete set of orthogonal 
vector wave functions to solve the structure's vector wave equation. In the two semi-infinite bottom and top 
media, those vector wave functions coincide with the plane-wave basis vectors, including both propagating and 
evanescent components. The technique is used to obtain the complete angular spectrum of the field scattered 
by a nanohole in a metallic film under Gaussian illumination. We also show how the obtained formalism can 
easily be extended to spherically and cylindrically multilayered media. In those cases, the expansion coefficients 
would multiply the spherical and cylindrical vector wave functions. 



I. INTRODUCTION 

With the last decade's revolution in microscopy and litho- 
graphic techniques, the field of nanophotonics has experi- 
enced immense growth. This expansion has prompted the de- 
velopment of new numerical and analytical techniques for the 
electromagnetic study of nanostructures. Some of the most 
important theoretical methods in nanophotonics are reviewed 
in 0] chap. 15]. The theoretical understanding of the interac- 
tion of light and matter at the nanoscale has allowed new dis- 
coveries in the fields of plasmonics, optical microscopy, and 
Raman scattering |(2]3- A particular problem in nanopho- 
tonics, which has particularly benefited from the new theoret- 
ical techniques, has been the extraordinary optical transmis- 
sion of light discovered by Ebbessen et al 0. This effect 
appears when we scatter light off back-illuminated subwave- 
length hole arrays in metallic thin films and has found appli- 
cations in sensing and lens design ((6] |7)). Since then, sev- 
eral authors have analyzed the transmission of light through 
isolated nanoholes in metallic thin films and the collective be- 
havior of arrays of nanoholes (see the extensive review in ]8) 
and references therein). 

In this work we present a technique which allows the plane- 
wave decomposition of the field scattered by an arbitrary per- 
turbation present in a planar multilayered structure under gen- 
eral external excitation to be obtained. The technique is ap- 
plied to the problem of a single back-illuminated nanohole in a 
metallic film. We show the exact angular spectrum representa- 
tion of the scattered electric field containing both propagating 
and evanescent components. To the best of our knowledge, 
this detailed field decomposition has not yet been studied for 
these nanostructures, even though the problem has received 
a significant amount of attention in recent years (for exam- 
ple, Refs. [9-13]). Understanding the spatial dependence of 
the scattered fields is necessary to efficiently collect the light 
transmitted through the nanoapertures, but more importantly, 
it can provide important information about the electromag- 
netic response of the nanostructure, which can then be used 



in applications like photovoltaic design [14], optical trapping 
ifTBll . and nano-optical tweezers fl6l . 

We finish the paper by outlining how the technique can 
be extended to spherically and cylindrically multilayered sys- 
tems to obtain scattered field decompositions in spherical and 
cylindrical vector wave functions. 

II. DERIVATION OF THE METHOD 

A monochromatic exp(— iwt) time dependence of the field 
has been assumed in this analysis. When faced with an elec- 
tromagnetic scattering problem, the dyadic Green's tensor 
(GT) technique can be used to solve it [17], offering some 
advantages over other tools like versatility and ease of cal- 
culation. Let us consider a fairly general scattering problem 
where a nonmagnetic base system, characterized by a space- 
dependent relative dielectric constant £{,(r), is illuminated by 
a general excitation Eo(r). Now, in some region of space 
V, we introduce a perturbation, i.e., the points r' £ V have 
a different dielectric constant with respect to the unperturbed 
system: 

e(r') = £ & (r') + Ae(r'). (1) 

Martin and Piller showed in ifTTl that the total electric field 
E(r) in the modified system can be expressed as: 

E = E + / fcg Ae(r') G (r, r')E(r')dr', (2) 
Jv 

where 

• Eo (r) is the electric field due to the external excitation 
that would be present in r if the perturbation did not 
exist, 

• fc 2 is the square of the vacuum wavenumber, 

• G (r,r') is the GT of the base (unperturbed) system, 
and 
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FIG. 1. (Color online) Geometry of the unperturbed system where 
the Green tensor is calculated. We consider a nonmagnetic planar 
multilayered structure, which is homogeneous in the (x, y) plane. 
The positions of the interfaces are given by di and the relative di- 
electric constant of each slab is given by ei . 



• E(r') is the total electric field present inside the pertur- 
bation. 

Clearly, the integral of the right-hand side of Eq. |2]) corre- 
sponds to the field scattered by the perturbation. One of the 
difficulties in applying this method is that obtaining the GT 
for an arbitrary base system is a daunting task. Fortunately, 



when the system of interest has strong symmetries like strat- 
ification, its GT can be calculated. References ifTSl . |[T9l and 
EOl provide the GT for planarly, spherically and cylindrically 
multilayered structures, respectively. In an excellent contribu- 
tion, Tan and Tan [21 1 developed a unified formulation for the 
three mentioned geometries which is also valid for biisotropic 
media. In all cases, the construction of the GT is based upon 
its expansion in the kind of vector wave functions appropriate 
for each geometry: rectangular, spherical or cylindrical vec- 
tor wave functions ETl sec. 4]. In this expansion, the GT is 
expressed as a sum of products of functions. Each term in the 
sum consists of the product of two functions, one which de- 
pends on r and the other on r'. This separate dependence is 
the key property exploited by our technique. 



A. Green's tensor for a Planar Multilayered Structure 

In this section we obtain a modal decomposition of the GT 
for a planar multilayered structure. The dyadic GT of a sys- 
tem represents the system's electromagnetic response to an 
infinitesimal excitation and allows computation of the elec- 
tric field induced in the system by an arbitrary distribution 
of electromagnetic sources [1, chap 2.10]. Let us examine 
the planar multilayered structure depicted in Fig. [T] Each 
layer I £ 1 . . . L has a different relative dielectric constant 
£; and relative magnetic permeability equal to one. The pla- 
nar boundaries between layers are perpendicularly oriented to 
the z direction and they are located at z coordinates di, I £ 
1...L—1. The layers extend to infinity in the x and y dimen- 
sions. 

Using IfTSl or fl2~Tl sec. 6], it is possible to calculate the 
effect of an electric dipole point source located at point r' = 
[x f , y', z'\ in a given layer I' on a target point r = [x, y, z] in 
layer /. Then we can arrive at the spectral decomposition of 
the structure's GT : 



S(ry) = -^(r-0 



+oo /*+00 



M„± 



'U,L(Mt,(r) ® >')) + N a% . k (N± (r) ® N* >')) 



(3) 



M feU y ( r ) = i(exp(±ik l z z) + s Rf; Kky exp(=fik l z z))s likaiky exp(i(k x x + k y y)) 
N feU„( r ) = ( e M±ikiz)pf;^ ky + p R^ ky exp(Tikiz)pf^ k J exp(i(k x x + k y y)) , 



(4) 
(5) 



where a ® b is the outer product of the three-dimensional 
vectors a and b; k\ = fc 0ll /£z are the wave numbers for 
I = 1...L; [k x ,k y ] are the transversal components of 
the wavevector, which are assumed to take only real values 
throughout the paper; v af; l k k for v — {M, N} are com- 
plex scalars which depend on the transversal wavevector, the 
source and target layer indexes, and the type of vector wave 



function (M or N); the rectangular vector wave functions 
k (r) and k (r) are the mode functions; | f2~2"l and for 
the ± selection, the signs on top are to be taken when z > z' 
and the bottom signs when z < z' . 

In Eqs. (4 1 and (5 1, k\ = J}tf - {k x + k%) where the 

branch with Im{/c^} > is chosen, and the definitions of the 
polarization vectors as a function of the transversal wavevec- 
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tor and the layer are, in the Cartesian basis, 

s _ {ky it. - k x y) 

U2 i 1.2 



kl + k% h 



(6) 



(7) 



These expressions coincide with the definitions of the S and P 
polarization vectors of a plane wave with transversal wavevec- 
tor [k x , k y ] propagating in a multilayered structure 11231 . 

Finally, we give the expressions for the expansion coeffi- 
cients: 



V(X ui k k 



(8) 



for v — {M, N}, or equivalently v = {s,p}. 



The v Rf k k and v Tfj k k quantities are sometimes re- 



ferred to as generalized Fresnel reflection coefficients and 
generalized transmission coefficients (see E4l chap. 3.3]). 
Recursive methods for computing v Rf k k and "T,/, k k are 
detailed in |T8), ED, ED, and E2 chap. 3.3]. 

Similarly to ll25l . we can establish orthogonality relation- 
ships between the vector wave functions. Given the following 
scalar product between two vectorial fields: 



(V|W> 



4-oo p-\-oo p-\-oo 



drV(r)W(r), (9) 



where V(r) is the Hermitian conjugate of V(r), and it is 
easy to verify that (M^ fe JM+ ?y ) = unless [k x ,k y ] = 



[qx,q y ]; (N+ fc JN+ gy ) = unless [k x ,k y ] = [q x ,q y ]; and 
<M+ tky |N+ t J = for all [k x , k y ] and [q x ,q y ]. 

The same orthogonality relationships apply for the minus 
superscripted vector wave functions. 

From now on, when applying ^ to the GT between a target 
point r = [x, y, z] and a set of source points r' = [x' , y' , z'] € 
Vj we will make the restriction 



either z > z' V(z, z') or z < z' V(z, z'). 



(10) 



As we will show shortly, this allows us to write our expres- 
sions in a very compact form. Also, for most scattering 
problems, the positions of electromagnetic sources and target 



points satisfy ( 10 1. In the next section this restriction will be 
applied to the positions of scatterers and target points as well. 



Taking ( 10 1 into account, we simplify the notation defining: 



J k x k. 



(r) = M fc , (r), e? fc »=N h u (r) 



f fe % (r') = M_ kx _ kv (r') , f* (r') = N_^(r'), 



where we drop the ± sign, which clutters the expressions 
and can easily be recovered with respect to ([3]) depending on 
sign(z — z'). We now re-write ^ with the newly defined 
vector fields: 



G(r,rO = -^(r-rO 



(11) 



E 



dk 



v<*in,Kk*l x k y {*)®tLky)- 



When the problem setup complies with restriction (flOj, we 
can use Eq. ( 1 1 1 to compute the field in a point r of the unper- 



turbed base system due to a source distribution J(r) present 
in volume Vj: 



E (r) 



dr' G (r,r') • J(r') 



/ d *'[- E / dk * / dk y a]>, LKky el ky (r) ® f^r')] ■ J(r') 

JVj K V v=sp J-oo J~oo 

^J(r) + ]T / dkj dk y et xky (T)( dr'a^ kxky fl ky (v')-3(r% 

l' v=s.pJ °° J— co JVj 



z 65 z 

~k 



z Q9 z 



v—s,p 



J(r)+J2 dk m dk y et xky (t)^ %Kky . 



(12) 



The first equality follows from using the GT dyadic G (r, r') 
to solve the inhomogeneous vector wave equation of the mul- 
tilayered structure. As indicated in [26 chap. 13.1], the GT 
can be used to solve such a system within a boundary surface 
S where the field must meet the required boundary conditions. 
The first line of Eq. ( |T2] > is obtained for S tending to infinity 



and enforcing that the field is zero at surface S. 

Outside volume Vj, where J(r) vanishes, the field is a 
weighted sum of the basis functions e kxky (r). The complex 
weight is given by: 

fe,= / dr'a^ kxk X ky (r').J(r'), (13) 
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The operation a ■ b represents either the product of a matrix 
and a vector or the inner product of two vectors. From Eq. 



( 12 1 we observe that, if we exclude the source points, the set 



k k (r) is a complete basis for the 



of vector wave functions et 
field distributions in the base structure. Since we have al- 
ready established their orthogonality, we conclude that, ex- 
cept at the source points, the set of vector wave functions 
e k k ( r ) i s a complete orthogonal basis for the field distribu- 
tions in the structure. Recall that depending on which side of 
restriction jfojl is met, ej^(r) = {M+ ky (r), N+ ky (r)} 

ore^^(r)={M fcsifev (r) I N^^(r)}. 

We can now comment on the convenience of using re- 
striction ( [Tol l. It i s c l ear that without imposing it we 
would get a spectral expansion with terms in both ba- 
^ {M+, fcB (r),N+ ky (r)} and {M" ky (r), (r)}, 
which would impair the clarity of the results and would be 
cumbersome to use. 



B. Dielectric Perturbations 

Now we are ready to get to the core of our technique. Let 
us go back to Eq. (|2|i describing the scattering problem when 
a perturbation with volume V and offset dielectric constant 
Ae(r) is introduced in the multilayered structure. The solu- 
tion of that equation is found in two steps: 

1. By discretizing the perturbation volume V in T small 
volumes with centers at r' = rt, and setting r = 
r t V t e {1 ... T}, a self-consistent 3T x 3T system 
of equations arises for the field inside the perturbation 
E(r t ). 

2. Once E(r t ) is available, the scattered field E sc (r) [in- 
tegral term of Eq. ^] and the total field E(r) can be 
computed for r £ V using E(r') as the generator in 
Eq. The volume integral becomes a finite sum over 
each discrete volume element. We keep the continuous 
notation for convenience. 

The first step is by no means trivial: see the introduction in 
E71 for a summary of the encountered difficulties. There are 
several possibilities to compute E(r') (UJ chap. 15]), but all 
of them typically have problems of convergence when the con- 
trast between the scatterer and the medium is high. In this pa- 
per we are interested in the expansion of the scattered field, so 
we do not discuss the subtleties of calculating the internal field 
E(r') for r' g V. For all purposes of this paper the internal 
field can be considered as an input to our problem. Just for the 
sake of completeness we mention that we calculated this field 
by using the self-consistent method given in 11281 . but note that 
our technique is totally independent of the particular method 
employed to calculate E(r'). 

We now focus on the integral part of the right hand side 
of Eq. (|2| when r V. As mentioned, this is the term that 
corresponds to the field scattered by the perturbation. Noting 
that e k k (r) does not depend on the integration variable r' 

and using the decomposition of G for r ^ r' in 



and the associativity of the vector outer and inner products, we 
arrive at our main result. By expressing the volume integral in 
((2) as 



, s; , -oo J-co 

we have hence arrived at a decomposition for E sc (r) in the 
vector basis modal functions of the unperturbed structure 

eLk» V r £ V. 



A restriction of the type ( 10 1 has been 
assumed in ( |14) . This time it applies to the positions of the 
scattering volume V with respect to the the target point(s) r. 
This assumption is again true in many scenarios of interest. 

We now turn our attention to the illuminating field Eo(r). 
Assuming that the target points we are interested in do not 
contain electromagnetic sources, Eq. (12i means that E (r) 
can be expanded in the e k k (r) basis by some complex 
scalars /3V,, k k . For the illuminating field, the source layer 
is the one where the active light source is. We define it to be 
the first one V = 1: 



Eo(r) = J2 



+ 00 



dk y f3l ukrk y k u(r). (16) 



The entire system of equations is linear in the incident field 
Eo(r). This can be easily checked by formally rewriting Eq. 
Q inside the perturbation volume V as an operator equation, 
i.e. E (r') = £E(r'). The integral operator C is linear. Thus, 
upon inversion the system remains linear, which means that 
the field inside the perturbation depends linearly on the inci- 
dent field Eo(r). It follows that by computing the decomposi- 
tion of the field scattered by the perturbation when illuminated 
by each of the e k k (r) basis functions, we can obtain the de- 
composition of E sc (r) due to an arbitrary excitation Eo (r) by 
mere linear superposition. 

Note that in the top and bottom semi-infinite layers, and 
because there are no reflections from ±oo, v Ri k k = 

v R~ Lk k — V [v,k x ,k y ]. This means that in those lay- 
ers e\ k k (r) and e v L k k (r) are just plane-waves [see Q 

and (|5j]. Note that both propagating (Im{fc] ,L } = 0) and 
evanescent (Im{fc^ L } ^ 0) plane waves are present in the 
expansion, giving us all the information contained in E sc (r). 







On those two layers, the type of expansion in equations ( 14 1 
and ( ff6| is commonly known as the angular representation of 
a field: (TJ chap 2.12] and (29] chap 3.2]. 

Let us now summarize the importance of results ( ff4p6 i. 
We formally start with an electromagnetic source J(r), away 
from the perturbation, which induces the field E (r) in the 
structure. This field is a weighted sum of the e]J k (r) func- 
tions. When the field E (r) encounters the perturbation 
Ae(r') in r' G V it interacts with it and, as a result, a scattered 
field E sc (r) is added to the original Eo(r). This new field is 
again a weighted sum of the ejr! k (r) functions when r ^ V. 
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In less abstract terms, assume for instance that an S (or P) 
linearly polarized plane wave illuminates an unperturbed pla- 
nar multilayered structure. The [k x ,k y ] components of the 
plane wave and its polarization will correspond to one of the 
e k k (r), which will exist in the structure without ever chang- 
ing [k x , k y ] or v. On the other hand, if there is a perturbation 
in the structure, the original plane wave will give rise to an 
infinite number of new [v,k x ,k y ] components with weights 



given by ( 15 1 



III. APPLICATION OF THE METHOD TO A NANOHOLE 

We envision many possible uses for this technique. In 
this section we will apply it to the well known problem of 
a nanohole in a metallic film illuminated by a Gaussian beam. 
In this way, we will show how to use this method in a typical 
problem and, at the same time, we will calculate the shape of 
the scattered field from a nanohole, which may be of impor- 
tance for certain applications. 
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FIG. 2. (Color online) Setup for nanohole scattering analysis. 
The excitation is a circularly polarized Gaussian beam propagating 
through the air side and impinging on the gold layer with zero in- 
cidence angle. The excitation wavelength is 950 nm and its Gaus- 
sian beam waist 1.5 ^im. The sample consists of a 400-nm-diameter 
nanohole made on a 140-nm layer of gold deposited on a 200- 
nm glass substrate. In the figure e\ = £4 = 1, £2 = £Au = 
-37.2 + 2.5i, and £3 = e g i ass = 2.25. Also di = 0, d 2 = -140 
nm, and ^3 = —340 nm. 

Fig. [2] shows a sketch of the considered system. The angu- 
lar spectrum decomposition is calculated at the air layer oppo- 
site of where the excitation comes from. The laser layer is the 
first one, the perturbation is contained in the gold layer I' = 2 
and the target layer is I = 4. Here is a step by step description 
of how the method is applied. 

To represent the Gaussian illumination, a set of basis func- 
tions e k k (r) is selected which needs to be sufficiently dense 
in [k x , k y ]. We call it i\ Since all its components are propa- 
gating in vacuum, the set can be restricted to functions whose 
tranversal wavevector lies inside the circle k$ = k x + ky. 
We can then use the plane wave expansion of a Gaussian 



beam from [30] in order to compute the expansion coefficients 
P\\ k k m equation ( 16 1. The relationship to the expansion 



coefficients in other layers is: 



PL 



(17) 



The interior of the nanohole is discretized at points rt. We 
assume that the field E^ fc (rt) present inside the perturbation 
when the external excitation is equal to each of the modes 
el t (r) G I is available. In our case, we have used a self- 
consistent dipole method to obtain it. 

We will define now the output basis set O of vectors 
e k x k y ( r ) over which the angular spectrum will be computed. 
The basis vectors in I and in O can be different, in order 
to properly describe the different properties of the incoming 
beam and the scattered field. In particular, if we are interested 
in the evanescent as well as the propagating components of 
the scattered field, the set O has to be chosen accordingly. 

Let us now select one of the E£ k (rt). Using ( 15 1 with the 
volume integral approximated by a discrete sum over the mesh 
points r t we can obtain the expansion coefficients 7^ k for 
each element of O. At the target layer, this set of coefficients 
is the angular spectrum of the scattered field if the excitation 
was only the mode e k k (r) £ I. This is repeated for every 
element in I and all the partial contributions are added with 
weights (3^ k k to obtain the final result. 

Should the ouput set O be a dense enough sampling of the 
[k x , k y ] space, a trivial integral in the momentum space would 
recover the value of the total electrical field E sc (r) at any 
point r of the structure meeting ( 10 1 by just using ( fl~4j>. 



Let us now analyze the results. Figs. 3(a) and |3(b) 



are 

three dimensional plots of It]^ | and \~/ k k | respectively. 
The clear resulting cylindrical symmetry of the scattered com- 
ponents is expected since structure, nanohole and illuminating 
beam have that symmetry. Note though, that none of the cal- 
culations in our method make use of the nanohole or beam 
symmetries: the method is independent of the shapes of the 
perturbation(s) and the illuminating beam, and its derivation 
does not contain any approximation. The strong presence of 
z-evanescent components illustrates the mechanism by which, 
through the volume integrals ( |15} , the perturbation generates 
components with transversal wavevectors absent in the origi- 
nal external excitation. When the volume of the perturbation 
or its dielectric contrast tend to zero, so do the amplitudes 
of these new components. Fig. |3(c)| is the three dimensional 
plot of \f3f 4 k k |, the angular spectrum representation of the 
S component of the direct field distribution at the target layer 
as written in ( 16 1 and (jlTJ. The corresponding plot |/3^ 4 k k \ 
is not includeasince it is practically identical. The direct field 
distribution is similar to the input Gaussian beam. It is mostly 
contained in the central region of the spectrum and decays 
rapidly with increasing normalized transversal wave vector 

k p = ^Jk x + ky/k . This similarity is due to the fact that the 

generalized transmission coefficients v Tr A , , do not have 
big variations across the central [k x , k y ] region for this partic- 
ular system. It is also apparent from Fig. [3]that there is much 
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FIG. 3. (Color online) Logarithmic plots of the absolute value of the expansion coefficients Jk^ky (a) and 7J! k (b), corresponding to the S 
and P polarization components of the angular spectrum of the scattered field and fti^k^k ( C X corresponding to the S polarization component 
of the angular spectrum of the direct field. The amplitudes are normalized to the total power (Po) of the input beam in layer I = 1, i.e. before 
going through the metal layer. The momentum axes are normalized to the vacuum wavenumber. 



more scattered light than direct light. Indeed, a simple calcu- 
lation reveals that the ratio of propagating power is about 250 
times in favor of the scattered field. 

Due to the aforementioned cylindrical symmetry, radial 
cuts of the figures contain the same information as the three- 



dimensional (3D) plots. Figures 4(a) and 4(b) are radial cuts 



of the 3D angular spectrum coefficients as a function of k p . 
The peaks in Figs. 4(a) and 4(b) correspond to normalized 
transversal wavevectors k p of 1.03 and 1.38. Both are evanes- 
cent in vacuum and their transversal wavevectors can be re- 
lated to two resonant modes of the multilayered structure PP . 
The P resonance position approximately matches the expected 
wavenumber of a surface plasmon at the gold-glass interface 
for the four layer system. It is interesting to see that a reso- 
nant S mode also appears. Even though it has a significantly 
smaller amplitude than the P mode, the difference in transver- 
sal wavevector means that the S mode field decay in the z di- 
rection is much slower than that of the P mode. The existence 
of these S-polarized waves in multilayered systems and their 



role in the extraordinary optical transmission of light have al- 
ready been discussed in the literature ll32l . 



A. Other Field Functional 

In this section, we have used different functions of the sets 
of coefficients 7^ k and 7^ k to expose some features of 
the angular spectrum representation of the scattered field of a 
nanohole under Gaussian illumination. It is worth highlight- 
ing that all the information contained in E sc (r) is captured by 
those coefficients. Consequently, any functional of the field 
must be obtainable from them. Let us take for instance a com- 
mon functional used in far field scattering analysis: the differ- 
ential scattering cross section. Following closely lfL2l sec. V], 
we define it here as: 



da 



^|E S c(r) 



(18) 



7 



-5.5 



-6.5 



-7.5 



7 log 10(^^) 



(a) 




-5.5 



-6.5 



-7.5 



0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.. 



(b) 



-loglO(^) 



0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.: 

fen 



FIG. 4. Logarithmic radial cut of (a) the S polarization components, 17^^ | and (b) the P polarization components |7j£ k 
field, with respect to the normalized transversal wavevector k p = ^k% + k^/kg. 



of the scattered 



where point r is supposed to be in the far field and rj depends 
on |r| 2 and the power of the illuminating beam. The differen- 
tial scattering cross section is a function of the solid angle 57, 
and it is related to the far field electromagnetic power that the 
perturbation scatters through each solid angle. Alternatively, 
it is a function of two linear angles, polar and azimuthal. In 
the far field these angles can be computed as 9 = arcsin(/c p ) 

and cf> — arctan (j^-^J respectively. The stationary phase ap- 
proximation l29l chap. 3.3] can be used to show that, at a 
point r = [a;, y, z] in the far field: 



|E sc (r)| 2 = C|fc z | 2 (| 7 ? k f + hl k f) 



(19) 



— y 



and C is constant 



with 4±1 = and h 

in the variables of interest. 

In our case, due to the cylindrical symmetry of the results, 
( |T~8| > will only depend on one of the angles, the polar angle 
= arcsin(fcp). To include only propagating components, it 
suffices to restrict k p < 1. Finally, to obtain 4|, rather than 
we make use of the fact that in our coordinate system 
dfl = d(pd0 cos(<9). 

For reference purposes, we provide the differential scat- 
tering cross section of the nanohole in Fig. [5] together with 
that of an infinitesimal electric dipole. Since our illuminat- 
ing beam has right circular polarization, it is appropriate to 
choose the orientation of the dipole to match it. The dipole 
is oriented along the x — iy direction. Its angular spectrum 
can be computed following Q] chap. 2.12.1], and the previ- 
ously outlined steps can be used again to obtain its differential 
scattering cross section. 

It can be noticed that the behavior of the scattered light from 
the nanohole is rather different than that of a single dipole. 
This is expected due to two main reasons: (a) on the one hand, 
in our case the size of the nanohole is only a fraction of the 




FIG. 5. (Color online) Differential scattering cross sections of the 
nanohole (cross-solid black line) and a x — iy oriented dipole (blue 



solid line). Both lines are normalized to their values at k p 



0. 



wavelength of radiation, meaning that we expect our radia- 
tion to be multipolar, and (b) on the other hand, the coupling 
to the surfaces modes greatly affects the radiation from the 
nanohole, as can be seen in Figs. |3]and|4] 



IV. CONCLUSION 

In summary, we have presented a method to obtain a plane 
wave expansion of the field scattered by a perturbation present 
in a planar multilayered structure under general illumination. 
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Once the field E(r') inside the perturbation is known, this 
method provides a very compact expansion of the scattered 
field, which allows for a straightforward analysis. 

The technique allows access to all the information con- 
tained in the scattered field by providing the amplitude and 
phase of the complex coefficients multiplying each plane wave 
for both orthogonal polarizations. This complete information 
can be easily manipulated to obtain other more common field 
dependent functions like the far field radiation diagrams and 
differential scattering cross sections, relevant when studying 
the far field properties of the light scattered by perturbations 
in a planar multilayered system. With respect to near-field 
applications like surface-enhanced Raman scattering or near- 
field microscopy, the technique should prove useful since it 
provides all the information about the evanescent components 
of the field, which are critical for these applications. 

In this article, the method has been applied to the problem 
of a nanohole in a thin gold film illuminated by a realistic 
Gaussian beam. The exact angular spectrum of its scattered 
field has been reported. Il33l . The diagrams present two dis- 
tinct peaks in the evanescent region, one for each polariza- 
tion. Both peaks can be traced back to resonant modes of the 
multilayered structure. 

Our technique can be easily extended to problems with dif- 
ferent geometries. From Eq. ([3]) on, our analysis is specific to 
planar geometries stratified in the z direction. Nevertheless, 



Ref. ||2T1 sec. 6] contains the Green's function decomposition 
of spherical and cylindrical multilayered structures as well, 
with f and p stratification directions. The steps in the pre- 
vious derivation can be repeated to obtain decompositions of 
scattered fields in cylindrical and spherical vector wave func- 
tions. In the case of spherical geometry the expressions will be 
formally identical to those obtained in Sec. [II] For the cylin- 
drical geometry, according to [20| and [21], there will gen- 
erally be cross terms in the expansion of the GT of the kind 
M(r) <g> N(r') and N(r) <g> M(r') (seen explicitly in 11201 ). 
This is not a blocking point. Since the separate dependence in 
r and r' is maintained, our derivation can still be carried out 
by introducing more terms besides e and / if necessary. 

We advance that this method can have an important impact 
in elucidating the different mechanisms of electromagnetic 
interaction between different features of nanostructures. The 
flexibility of extending this method to other geometries will 
prove very useful in a large number of applications. 
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